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I.  INTRODUCTION 


A  number  of  investigations  have  provided  evidence  that  simple, 
spherically  symmetric  source  representations  are  inadequate  to 
explain  surface  wave  generation  by  underground  explosions.  For 
example,  Murphy  (1977)  demonstrates  that  a  spherically  symmetric 
source  model  subject  to  the  depth/yield  scaling  proposed  by  Mueller 
and  Murphy  (1971)  is  consistent  with  a  wide  variety  of  near-regional 
and  teleseismic  body  wave  data;  he  shows,  however,  that  the  same 
model  is  inconsistent  with  long-period  surface  wave  observations. 
In  particular,  the  source  scaling  relationships  predict  that  surface 
wave  magnitude  will  increase  with  yield  (W)  less  rapidly  than  log  W, 
whereas  studies  by  Basham  and  Horner  (1973),  Springer  and  Hannon 
(1973)  and  others  consistently  give  yield  exponents  exceeding  1.0 
(about  1.1  to  1.3).  A  second  example  pertains  to  anomalies  in  the 
Rayleigh  waveform,  rather  than  amplitude.  Rygg  (1979)  demonstrates 
that  the  Rayleigh  waves  from  an  eastern  Kazakh  explosion  exhibit  a 
polarity  reversal  and  time  shift  relative  to  similar,  nearby 
explosions. 

In  both  these  cases,  the  authors  proposed  that  spall  closure 
was  responsible  for  the  anomalous  surface  wave'  observations.  To 
quote  Rygg,  "The  only  source  that  appears  capable  of  introducing 
both  phase  reversal  and  time  delay  is  the  phenomenon  of  spall 
closure."  Such  anomalous  surface  wave  observations  have  been  widely 
attributed  to  the  effect  of  spall  closure,  and  two  principal  reasons 
can  oe  cited.  Firstly,  spall  of  near-surface  material  is  a 
we 1 1 -documented  and  nearly  ubiquitous  component  of  contained 
underground  explisions  (for  example,  Eisler  and  Chilton,  1964; 
Viecelli,  1973;  Springer,  1974).  Secondly,  Viecelli  has  purportedly 
shown  theoretically  that  spall  closure  can  be  a  major  source  of 
Rayleigh  wave  radiation  from  buried  explosions. 

If  spall  were  the  primary  source  of  explosion-induced  surface 
waves,  or  an  important  secondary  source,  this  fact  would  have 
important  consequences  for  the  seismic  estimation  of  explosion 
yield.  Shots  with  anomalously  large  spall,  say  due  to  underburial. 
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would  in  turn  have  anomalously  large  surface  wave  excitation.  In 
order  to  reliably  exploit  surface  wave  magnitudes  in  the  seismic 
determination  of  explosion  yield,  it  is  essential  to  understand  the 
potential  effects  of  spall  on  the  surface  wave  signal. 

The  purpose  of  this  paper  is  to  point  out  and  correct  a 
fundamental  flaw  in  previous  theoretical  work  on  the  subject  of 
surface  wave  excitation  by  surface  spall.  By  means  of  a  simple 
physical  argunent,  it  is  demonstrated  that  spall  is  incapable  of 
generating  a  significant  contribution  to  long-period  (greater  than 
about  ten  seconds)  surface  waves.  It  is  shown  that  previous 
conclusions  to  the  contrary  have  been  based  on  source  descriptions 
which  do  not  conserve  momentum.  An  equivalent  elastic  point-source 
representation  is  derived  for  spall  which  corrects  this  momentum 
imbalance.  This  simple  source  model  predicts  that  spall  can'  enhance 
surface  wave  excitation  only  for  relatively  short-period  waves  —  on 
the  order  of  a  few  seconds  period  —  with  negligible  effect  on 
20-second  radiation. 

We  believe  that  the  argument  against  long-period  surface  wave 
enhancement  by  spall  is  conclusive,  since  it  rests  only  on  the  very 
fundamental  principle  of  conservation  of  momentum.  In  particular,  * 
the  conclusion  would  seem  to  be  independent  of  assumptions  about  the 
nonlinear  rheology  of  the  source  region.  There  remains  the  j 

possibility,  however,  that  the  occurrence  of  spall  might  indirectly 
modify  the  surface  wave  signal  by  altering  the  effective  reduced 
displacement  potential  (RDP),  or  monopole  component,  of  the 

explosion.  To  test  the  generality  of  our  conclusion,  therefore,  we 
compare  the  simple  model  with  predictions  from  a  nonlinear  numerical 
simulation  of  a  buried  explosion. 
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II.  REVIEW  OF  THE  PROBLEM 


Spalling  is  defined  (Eisler  and  Chilton,  1964)  as  the  parting 
of  •  near-surface  layers,  which  were  originally  in  contact,  in 
response  to  the  stress  waves  produced  by  a  contained  underground 
explosion.  The  explosion-generated  stresses  become  tensile  upon 
reflection  at  the  free  surface,  and  spall  results  when  the  tensile 
strength  of  the  near-surface  geologic  material  is  exceeded.  The 
detacned  layer  then  decelerates  under  the  force  of  gravity,  and 
eventually  slaps  down. 

It  is  the  slapdown  phase  that  has  been  the  focus  of  attention 
as  a  potential  generator  of  surface  waves.  Near-field  acceleration 
pulses,  radiated  upon  slapdown  ("spall  closure")  can  be  identified 
on  free-surface  accelerometer  recordings  (Eisler  and  Chilton,  1964; 
Rimer,  et  al_. ,  1979),  from  which  the  shape  and  duration  of  the  spall 
gap  can  be  estimated.  Eisler  and  Chilton  have  analyzed  travel  time 
data  from  such  recordings,  from  which  they  deduced  that  the  shape  of 
the  spall  gap  is  roughly  conical,  with  sides  concave  upward.  For 
large  explosions,  ballistic  periods  as  long  as  two  seconds  are 
inferred  (Viecelli,  1973). 

Viecelli  (1973)  studied  surface  accelerometer  records  from 
nuclear  explosions  in  order  to  estimate  the  mass  and  momentum  of  the 
spall  region.  He  determined  the  following  approximate  empirical 
relationships  between  the  explosion  yield  (W),  and  the  spall  mass 
(m$)  and  momentum  (I$): 

ms  *  1.6  x  109  W  Kg  (1) 

Is  *  4.6  x  109  W  Nt-sec  ,  (2) 

where  W  is  in  kilotons.  These  estimates  formed  the  basis  for  his 
analysis  of  the  surface  wave  contribution  from  spall  closure. 
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To  estimate  the  surface  wave  excitation,  Viecelli  approximated 
the  slapdown  phenomenon  as  a  vertical  load  F  applied  at  the  free 
surface.  The  time  history  was  assumed  to  be 

F(t)  -  V0Mt2  +t2)  , 

where  tq  is  a  characteristic  time  representing  the  duration  of  the 
spall  impact  (about  twice  the  spall  thickness  divided  by  the  P-wave 
velocity).  Actually,  since  Xg  is  short,  on  the  order  of  a  few 
tenths  of  a  second  at  most,  and  the  surface  wave  periods  of  interest 
are  many  times  longer,  this  time  history  can  be  viewed  as 

essentially  a  delta  function  with  area  I$.  Using  the  analytic 
solution  for  the  surface  waves  from  a  vertical  load  on  a  uniform 
elastic  halfspace,  and  the  empirical  estimate  for  I$,  Equation 
(2),  it  was  found  that  this  impulsive-load  model  predicts  surface 
wave  amplitudes  comparable  to  those  observed  for  underground 
explosions.  Even  without  adding  the  direct  contribution  from  the 
explosive  monopole,  it  was  argued,  the  spall  momentum  is  sufficient 
to  produce  surface  waves  of  the  observed  amplitude.  Thus,  spall 
would  appear,  on  this  basis,  to  be  a  potential  contributor  to 
anomalous  surface  wave  amplitude  and  phase  observations. 

A  different  theoretical  approach  to  estimating  the  spall 

effect  on  surface  waves  was  taken  by  Harkrider  and  Bache  (1979).  We 
can  refer  to  this  approach  as  the  wave-suppression  method.  Analytic 
methods  for  computing  synthetic  surface  wave  seismograms  for  buried 
sources  in  layered  elastic  media  represent  the  source  by  means  of  a 
stress-displacement  discontinuity  at  the  source  depth  (for  example, 
Harkrider,  1964).  This  representation  is  equivalent  to  specifying 
the  up-  and  downgoing  waves  emitted  at  the  source.  The  surface 
waves  associated  with  each  of  these  components  can  be  calculated 
separately,  as  described  by  Harkrider  and  Bache. 

In  the  wave-suppression  approach,  the  effect  of  spall  is 

modeled  by  partial  or  complete  suppression  of  the  source 

contribution  from  upgoing  waves  generated  by  the  explosion.  This 


4 


SYSTEMS.  SCIENCE  A  NO  SOETW 


reduction  of  the  upgoing  source  wavefield  is  intended  to  simulate, 
in  an  approximate  way,  the  loss  of  energy  from  the  upgoing  wave  due 
to  spallation  or  other  near-surface  nonlinear  phenomena.  Bache, 
Goupillaud  and  Mason  (1977)  used  this  technique  to  model  spall  for 
explosions  in  NTS  source  materials.  They  found  that  complete 
suppression  of  the  upgoing  waves  reduced  the  predicted  M$ 
(measured  at  about  20  seconds  period)  by  about  0.4  units  for 
explosions  in  granite.  For  explosions  in  tuff,  the  effect  was  to 
enhance  M$  by  about  0.1  to  0.2  units.  Harkrider  (1981)  points  out 
that,  for  an  explosion  in  a  soft  material  overlying  a  hard  material, 
suppression  of  the  upgoing  wave  can  reverse  the  polarity  of  the 
Rayleigh  wave. 

Both  these  theoretical  approaches,  the  impulsive-load  model 
and  the  wave-suppression  model,  predict  that  spall  can  substantially 

modify  the  surface  waves  generated  by  a  contained  explosion.  The 

actual  physical  processes  involved  in  spall  are  complex,  and  the 

predictions  of  the  simple  elastic  theories  are  fairly  sensitive  to 
assumptions  which  are  difficult  to  test.  For  example,  there  is  no 
physical  basis  for  determining  what  fraction  of  the  upgoing  waves  to 
use  in  the  wave-suppression  model. 

For  this  reason,  numerical  explosion  simulations  should  be  an 
important  theoretical  tool  for  evaluating  the  spall  effect  and 
testing  the  elastic  theories.  Viecelli  (1973)  compared  two  finite 
difference  simulations,  one  in  which  tension  failure,  and  therefore 
spall,  was  permitted,  and  a  second  in  which  only  elastic  behavior 

was  permitted.  Gravitation  was  included,  so  the  spall  model 

included  closure.  The  computed  Rayleigh  wave  amplitude,  observed  in 
the  finite  difference  grid  at  a  range  of  15  km,  was  2.7  times  larger 
for  the  spall  model  than  for  the  elastic  model,  apparently 

corroborating  the  results  of  the  elastic  theory.  An  important 
feature  of  the  result,  however,  is  evident  in  Viecelli 's  Figure  10: 

the  Rayleigh  wave  amplitude  has  been  measured  at  a  predominant 

period  of  about  three  seconds.  We  will  remark  on  this  further  in  a 
later  section. 
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Both  the  simple  elastic  theories  and  Viecelli's  numerical 
simulation  agree  in  predicting  a  potentially  important  spall  effect 
on  teleseismic  surface  waves.  Yet,  we  find  in  the  next  section  that 
the  elastic  theories  proposed  so  far  are  flawed,  and  their  results 
consequently  are  misleading.  The  results  of  Viecelli's  numerical 
simulation,  on  the  other  hand,  should  simply  be  reinterpreted.  In  a 
later  section,  we  analyze  a  new  simulation  in  terms  of  its 
teleseismic  surface  wave  radiation,  using  the  procedures  of  Bache, 
Day  and  Swanger  (1982).  These  results  will  be  found  to  compare 
favorably  with  a  corrected  elastic  theory. 
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III.  ANALYSIS  OF  THE  ELASTIC  MODELS 


In  the  last  section  we  reviewed  two  elastic  source  models 
which  have  been  used  to  approximate  the  effect  of  spall  on  surface 
waves.  Although  the  physical  process  of  spall  is  highly  nonlinear, 
it  is  a  valid  analytical  approach  to  approximate  it  by  elastic 
loading.  It  is  well  known,  of  course,  that  in  general,  elastic 
sources  exist  whose  radiation  is  exactly  equivalent  to  a  given 
nonlinear  seismic  source.  Our  critique  of  the  proposed  elastic 
models  centers  on  their  failure  to  conserve  momentum.  It  will  be 
shown  that  when  this  deficiency  is  corrected,  a  simple  elastic  model 
of  spall  predicts  negligible  surface  wave  excitation  at  periods  of 
teleseismic  interest. 

For  now,  we  avoid  tying  the  analysis  to  any  assumptions  about 
the  source-region  rheology  by  imagining  the  source  region,  V$,  to 
be  enclosed  by  a  surface,  E,  lying  wholly  outside  the  zone  of 
nonlinear  response  (Figure  1).  The  surface-wave  generation  can  be 
expressed  in  terms  of  the  elastic  properties  outside  E  and  the 
displacements  (ua)  and  surface  tractions  (ta)  on  E.  The 
displacements  and  tractions  on  E  do  depend  on  the  rheology  of  the 
source  region,  but  we  will  be  able  to  sidestep  this  problem 

temporari ly. 

Bache,  Day  and  Swanger  (1982)  give  an  expression  for  Rayleigh 
wave  spectral  displacement  u~?  in  the  low-frequency  limit,  for  a 
general  axisymmetric  explosion  source.  It  is  assumed  that  the  only 
body  force  acting  on  the  source  volume  is  gravity,  and  that  the 
source-generated  displacements  and  tractions  are  taken  relative  to 
their  initial,  static-equilibrium  values.  Then  we  have 


1 1m  U  R(r,z) 
a 


S?  (r,z)  [I 

3  I 


*  V 


where 


(3)  | 
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Figure  1.  Geometry  for  the  surface-integral  representation  of  the 
seismic  radiation  from  an  axisymmetric  explosion  source. 
The  receiver  point  is  (r,z),  and  a  typical  point  lying  on 
the  imaginary  surface  E  is  (r0,z0).  The  surface 
lies  wholly  outside  the  nonlinear  region  surrounding  the 
source. 
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The  notation  is  as  follows: 

a  *  cylindrical  component  (r  or  z) 

cR  »  phase  velocity  of  the  Rayleigh  wave  mode 

eQ  »  surface  ellipticity  of  the  mode 

x,u  *  Lame  constants  on  Z 
r0,Zg  »  cylindrical  coordinates  on  I 

r,z  »  cylindrical  coordinates  of  the  receiver  point 

n^  »  component  of  the  unit  normal  to  Z 


In  'S’j,  (r,z),  the  factors  w/wQ,  u*/wq,  and  AR  are  as  defined 
by  Harkrider  (1964):  the  normalized  vertical  and  horizontal 

eigenfunctions  and  the  source-independent  amplitude  response  of  the 
>  structure,  respectively;  and  are  Hankel  functions.  The 
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factor  <ifa  depends  only  on  the  layered  elastic  structure  outside  E, 
and  the  integrals  1^  and  I ^  depend  on  the  source-generated 
displacements  and  tractions  on  E.  At  times  long  compared  to  the 
source  duration,  the  displacements  and  tractions  reach  some  static 
values  on  E,  so  in  general  their  transforms  u  and  t  are 
proportional  to  u»  at  low  frequency.  Therefore,  ^(u)  is  of 
order  1: 

I2  “  0(1)  . 

The  1^  term,  on  the  other  hand,  is  simply  the  spectrum  of  the 
total  vertical  force  on  E (minus  its  initial  value  in  equilibrium 
with  gravity).  Conservation  of  momentum  demands  that  the  time 
integral  of  the  total  vertical  force  on  z  go  to  zero  at  times  long 
compared  to  the  source  duration;  this  is  true  irrespective  of  the 
nonlinearities  inside  E.  Therefore,  1^  goes  to  zero  at  least  as 
fast  as  at : 


1^  ~  0(<*ir)  ,  v  >  1  . 

The  source  duration  will  be  dominated  by  the  ballistic  period  of  the 
spall,  which  is  at  most  a  few  seconds.  So,  for  periods  of  interest 
for  teleseismic  surface  waves,  we  expect  I2  to  dominate,  being  of 
order  1,  and  1^  should  be  negligible. 

Now  the  unacceptability  of  the  impulsive-load  elastic  model  is 
apparent.  The  explosive  monopole  imparts  zero  net  momentum  to  the 
continuum,  and,  correctly,  makes  a  contribution  of  order  1  only  to 
Ig.  The  impulsive  surface  load,  however,  imparts  a  net  downward 
momentum  equal  to  Is,  the  spall  momentum.  Thus,  for  the 
impulsive-load  model,  we  have 

IL(0)  -  Is  ; 
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so  the  I ^  term  is  of  the  same  order  as  the  ^  term.  The  large 
spall  contribution  predicted  by  the  impulsive  load  model  is 
therefore  directly  attributable  to  its  inherent  momentum  imbalance, 
which  cannot  occur  for  any  physical  model. 

The  wave-suppression  model  is  similarly  flawed  by  failure  to 
conserve  momentum.  Lixe  the  impulsive-load  model,  the 
wave-suppression  model  imparts  a  net  downward  momentum,  again 
resulting  in  a  spurious  long-period  contribution  to  1^.  To  see 
the  momentum  imbalance,  we  examine  the  source  wavefield  radiated  by 
a  monopole  explosion  source.  The  source  strength  is  defined  by  a 
reduced  displacement  potential  'f(w),  such  that  the  radiated  spectral 
displacement  u-  is 


u(r,z,«)  -  -  7  fTl  e1  “  R/a  ? 


?  ?  1  /  ? 

where  a  is  the  P-wave  speed,  R  =  (r  +  z  )  ,  and  we  have 

taken  the  coordinate  origin  to  coincide  with  the  source  point. 
Expanding  the  term  in  brackets  via  the  Sommerfeld  integral  gives 


/oo 

dk  k  sgn(z)  JQ(kr)  z  - 


(kr)  ?  e  a 


where  _r  and  i  are  unit  vectors  in  the  radial  ana  vertical 
directions,  respectively,  and 


/  2,  2  ,2,1/2 
v  n  ( a  -  K  J 
a 


From  Equation  (4)  we  can  obtain  the  traction  vector  7  in  the  z  plane: 


u  a 

♦  sgn(z)  2 ;ik  J1(kr)r  eiva  I  2  •  . 

- 
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The  vertical  traction  is  continuous,  so  no  net  vertical  force 
acts,  and  momentum  is  conserved.  However,  suppressing  a  fraction  n 
of  the  upgoing  wave  requires  multiplying  the  exponential  factors  in 
Equations  (4)  and  (5)  by  [1  -  Q.5n  (1  -  sgn(z))]  (taking  z  positive 
downward).  The  result  is  to  introduce  a  discontinuity  in  the 
vertical  traction  at  z  *  0,  with  Fourier-Bessel  transform  «tz 
given  by 


a?z(k,»)  »  nu  (i  vj"1  (2k2  -  «2/e2)w(»)  .  (6) 

The  total  momentum  I  due  to  this  discontinuity  in  traction  is  just 
proportional  to  at  zero  frequency  and  zero  wavenumber: 

I  -  -  2»  «t2(0,0) 

(7) 

■  n  2i  pa  f  , 

OO* 

where  f^is  -iuf(O),  the  late-time  limit  of  't(t). 

Thus,  the  wave-suppression  model  does  not  conserve  momentum, 
and  the  momentum  imbalance  makes  a  spurious  low-frequency 
contribution  of  order  1  to  1^  of  Equation  (3).  The  net  momentum 
imparted  is  downward.  Using  Murphy's  (1977)  value  of  f*,  =>  4.2  x 
10  m^  for  a  100  KT  event  in  tuff,  the  magnitude  of  the  momentum 
imbalance  for  the  wave-suppression  model,  with  fractional 
suppression  n  as  1/4,  is  approximately  equal  to  Viecelli's  estimate 
of  the  slapdown  momentum. 


12 


srsrrMS.  science  a  no  softwarc 


IV.  AN  ELASTIC  SOURCE  MODEL  FOR  SPALL 


The  impulsive-load  and  wave-suppression  models  of  spall  are 

both  unphysical,  in  that  each  imputes  a  net  momentum  to  the  source. 
Furthermore,  long-period  surface-wave  predictions  are  seriously 
influenced  by  this  momentum  imbalance,  as  shown  by  Equation  (3). 
Here  we  consider  a  very  simple  model  of  spall  and  construct  its 

equivalent  elastic  loads.  It  will  be  seen  that  this  model  predicts 

2 

a  spall  contribution,  proportional  to  w  I  ,  which  is  certain  to 
be  negligible  at  20  second  periods. 

The  impulsive-load  model  reviewed  earlier  accounts  for  the 

inertial  forces  delivered  to  the  intact  continuum  by  the  spall  mass 
upon  closure.  Its  failure  to  conserve  momentum  is  due  to  neglect  of 
(1)  the  inertial  forces  imparted  during  spall  opening  and  (2)  the 

relaxation  of  gravitational  forces  during  the  ballistic  period  (we 
take  the  initial  equilibrium  configuration  as  the  reference  stress 
state).  If  the  spalled  material  is  severed  at  an  average  upward 
velocity  Vg,  and  if  the  period  over  which  the  mass  accelerates  to 
Vq  is  negligible  compared  to  the  periods  of  interest,  we  can 

approximate  the  inertial  force  of  spall  opening,  f^,  by  a  downward 
delta  function, 

fx-  I$  «(t) 

where  I$  is  the  spall  momentum  msVg.  Subsequent  to  the  spall 
opening,  the  continuum  experiences  a  relaxation  of  the  gravitational 
forces  exerted  by  the  spall  mass,  f2,  given  Dy 

f2«  -  m$g  [H( t )  -  H(t-Ts)]  , 

where  mg  is  the  spall  mass,  g  is  the  gravitational  acceleration, 
and  T$  is  the  average  ballistic  period. 
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Finally,  the  slapdown  rapidly  applies  inertial  forces  with  time 
integral  I$,  which  we  can  again  approximate  by  the  delta  function: 

f3  -  Is  4(t)  . 

These  equivalent  forces  are  sketched  in  Figure  2.  Obviously,  when 
all  three  terms  f^  +  *  f3  *  f$  are  considered,  the 

equivalent  forces  correctly  conserve  momentum.  Our  elastic  source 
model  for  spall  consists  of  applying  these  equivalent  forces  as 
point  loads  at  ground  zero,  superposed  on  the  explosion  monopole. 
We  neglect  the  finite  spatial  distribution  of  the  equivalent  forces, 
and,  as  already  stated,  we  are  ignoring  the  finite  duration  of  the 
inertial  forces  of  spall  opening  and  closing.  This  approximation 
will  be  adeouate  for  periods  substantially  greater  than  the  time  for 
a  wave  to  cross  the  spalled  region,  say  one  or  two  seconds. 

The  Rayleigh  wave  vertical  displacement  due  to  a  monopcle  at 
depth  Zq  with  reduced  displacement  potential  ?(»)  is  (Harkrider, 
1964) 


IV  2)  .  2„i„<,2  «  A  H t2,/s:  1  fi.  “*(*0>  ♦ 

u2(r,z)  •  A*  H0  ^Cr  j  ^Cr  *q  „0 


and  the  Rayleigh  wave  due  to  the  vertical  surface  load  7  (w)  is 


“z(r.z)  -  “  i  7S  H02)  r  ) 

Therefore,  the  presence  of  spall  modifies  the  monopole  Rayleigh  wave 
spectrum  by  the  factor 


+  II  (4,p  a2)’1 
? 


4*(zn)  *'(zn) 


X 


X 


-1 


(8) 
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fj  s-msg  [H(l)-H(t-T,)] 
fs»I,8(t-T,) 


fi  +  fp  -rf-! 


f  dt  (f  j 


Figure  2.  Equivalent  elastic  surface  load  for  an  idealized  spall 
model.  Correctly  accounting  for  both  inertial  and 
gravitational  forces  associated  with  the  spall  mass 
insures  zero  net  source  momentum. 
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where 


and  e  is  the  delay  between  detonation  and  spall  opening,  which  we 
will  approximate  as  zn/a.  The  fractional  contribution  of  the 
spall,  A-l,  is  therefore  proportional  to  w  at  low  frequency: 


A-l  ~  0(«2)  , 

(since  f  is  normally  of  order  w"*).  This  result  is  in  contrast  to 
the  impulsive-load  and  wave-suppression  models,  for  which  we  have 
seen  that  A-l  -  0(1). 

Figure  3a  shows  the  surface  wave  enhancement  factor.  A,  as  a 
function  of  frequency,  for  several  different  values  of  the  spall 
momentum  and  spall  mass.  The  f(w)  used  in  Equation  (8)  was  obtained 

from  a  numerical  simulation  of  a  spherically  symmetric,  61  KT 

4  3, 

explosion  in  granite  ('i'a0  »  1.6  x  10  m  ).  The  spectrum  of  the 
reduced  velocity  potential  is  shown  in  Figure  3b.  The  values  of 
I$  and  m$  as  given  in  Figure  3a  have  been  normalized  to 
Viecelli's  mean  empirical  estimates  for  61  KT  (Equations  1  and  2). 
In  no  event  does  spall  make  a  significant  contribution  to  the  20 
second  Rayleigh  wave.  The  maximum  enhancement  occurs  at  a  period  of 
about  1.4  Tg.  This  short-period  enhancement  might  be  important 
for  interpreting  regional  and  near-regional  surface-wave  data,  but 
will  have  no  effect  on  teleseismic  determinations  of  M$. 

Now  we  can  reexamine  Viecelli's  (1973)  numerical  simulation  of 
spall  in  light  of  the  foregoing  analytical  results.  Recall  that  the 
numerical  simulation  showed  Rayleigh  wave  amplitude  enhancement  of  a 
factor  of  2.7,  which  appeared  to  support  the  impulsive-load  elastic 
model.  However,  the  Rayleigh  wave  amplitude  was  measured  very  near 
the  source  (15  km),  where  the  dominant  period  of  the  computed  motion 
was  about  3  seconds.  From  Figure  3,  we  see  that  this  is  in 
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Case 


3.92  14 


1.96  5.6 
0.98  3.5 


Frequency  (Hz) 


Frequency  (Hz) 


Figure  3.  (A)  Theoretical  Rayleigh  wave  enhancement  factor  for  the 

equivalent  elastic  load  model  of  spall  described  in  the 
text  and  in  Figure  2.  The  and  Is  values  represent 
spall  mass  and  momentum,  divided  by  Viecelli's  (1973) 
mean  empirical  estimates  for  61  KT.  (B)  Reduced  velocity 
potential  for  a  61  KT  one-dimensional  explosion 
simulation,  which  was  used  in  Equation  8  to  construct  the 
Rayleigh  wave  enhancement  factors  shown  in  (A). 
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precisely  the  period  range  where  enhancement  is  predicted  by  the 
momentum-conserving  equivalent-force  model.  The  results 
Viecelli's  numerical  simulation  are  therefore  consistent  with  the 
equivalent-force  model,  and  need  not  be  interpreted  as  supporting 
the  impulsive-load  model. 
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V.  NONLINEAR  NUMERICAL  SIMULATION  OF  SPALL 


The  preceding  theoretical  argument  against  long-period  surface 
wave  enhancement  by  spall  appears  to  be  quite  conclusive,  since  it 
relies  only  on  the  very  fundamental  principle  of  conservation  of 
momentum.  In  particular,  the  conclusion  would  seem  to  be 
independent  of  assumptions  about  the  nonlinear  rheology  of  the 
source  region.  There  remains  the  possibility,  however,  that  the 
occurrence  of  spall  might  indirectly  modify  tne  long-period  surface 
wave  signal  by  altering  the  effective  reduced  displacement  potential 
of  the  explosion.  Such  a  phenomenon  was  envisioned  by  Aki,  Bouchon 
and  Reasenberg  (1974),  for  example,  who  suggested  that  spall 
slapdown  might  compact  the  nonlinear  source  volume,  leading  to  a 
reduction  of  the  long-period  level  of  the  RDP. 

To  test  the  generality  of  our  conclusion,  we  have  to  rely  on 
numerical  simulations  of  contained  explosions,  which  can  incorporate 
the  nonlinear  material  response  of  the  source  region.  Since  we  are 
interested  in  nonlinear  phenomena,  it  is  important  that  the 
simulations  preserve  relatively  short  wavelength  features  of  the 
near-source  disturbance,  even  though  our  ultimate  interest  is  in  the 
excitation  of  surface  waves  with  wavelengths  of  tens  or  hundreds  of 
kilometers.  To  model  the  spall  process,  the  source  simulation  must 
be  two-dimensional  (at  least),  so  as  to  include  the  free  surface, 
gravity,  and  the  depth  dependence  of  overburden  pressure;  these 
phenomena  are  absent  from  one-dimension  (spherically  symmetric) 
explosion  simulations.  Comoarison  of  two-dimensional  simulations 
with  one-dimensional  simulations  which  use  the  same  constitutive 
models  then  provides  a  means  of  quantifying  the  effects  of  spall  on 
the  seismic  radiation. 

The  finite  difference  method  is  suitable  for  treating  the 
two-dimensional,  nonlinear  problem,  but  it  is  not  feasible  to 
compute  the  wavefield  at  distances  beyond  a  kilometer  or  so  from  the 
source  by  this  method.  Provided  the  finite  difference  simulation 
extends  out  to  the  range  of  linear  material  response,  however,  the 
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teleseismic  signal  can  be  synthesized  from  the  close-range  finite 
difference  solution  using  Green's  function  methods.  This  approach 
was  taken  by  Bache,  Day  and  Swanger  (1982),  who  synthesized 
teleseismic  Rayleigh  waves  from  two-dimensional  explosion 
simulations  which  included  spall.  Comparing  synthetics  from  one- 
and  two-dimensional  explosion  simulations,  they  noted  that 
two-dimensional  source  phenomena  had  very  little  effect  on  either 
the  amplitude  or  waveform  of  the  long-period  teleseismic  Rayleigh 
wave,  but  they  observed  some  enhancement  of  very  short  period  (2  to 
4  seconds)  components.  This  result  is  in  excellent  agreement  with 
the  momentum-conserving  elastic  model  of  spall  introduced  in  the 
last  section.  However,  the  explosions  considered  by  Bache,  et  al. 
were  overburied  (depths  of  188  W*^3  meters,  compared  to  118  W^3 
meters  for  PILE0R1VER,  for  example);  so  we  still  cannot  exclude  the 
possibility  that  spall  from  shots  at  normal  containment  depth  (~  120 
W1'3  meters)  mignt  modify  the  long-period  surface-wave  signal. 

To  address  this  issue,  a  two-dimensional  ( ax i symmetric) 

simulation  of  a  buried  explosion  in  granite  was  performed  using  an 

\ 

explicit  finite  difference  method.  The  source-region  geologic 
structure  was  a  three-1 ayered  halfspace  representing  the  geology  of 
the  1966  PILEDRIVER  explosion  at  the  Nevada  Test  Site  (MTS).  The 
nonlinear  material  model  used  in  the  calculation  is  described  by 
Rimer,  et  al_.  (1979),  and  includes  tensile  failure,  shear  failure, 
irreversible  pore  crushup,  and  an  effective  stress  law.  The 
explosion  depth  and  yield  were  the  same  as  for  PILEDRIVER,  463 
meters  and  61  XT,  respectively  (i.e.,  depth  *  118  W^3  meters). 
The  initial  stress  field  was  taken  to  be  lithostatic. 

Figure  4  shows  the  computed  vertical  velocity  time  histories 
obtained  from  the  simulation  at  various  depths  directly  below  ground 
zero  and  at  various  ranges  along  the  free  surface.  These  are 
compared  with  recorded  velocities  for  the  PILEDRIVER  shot.  The 
spall  occurs  at  all  the  sites,  immediately  after  the  first  velocity 
peak,  as  evidenced  by  the  -l  g  slopes  of  the  velocity  waveforms. 
Spall  closure  can  be  identified  with  the  termination  of  the  -1  g 
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VELOCITY  (m/*) 


ZERO  RANGE  ZERO  DEPTH 


Figure  4.  Vertical  velocity  time-histories  (positive  up)  obtained 
from  the  two-dimensional  simulation,  compared  with 
recorded  velocities  for  PILEDRIVER.  The  data  are  from 
Hoffman  and  Sauer  (1969).  Occurrence  of  soall  is 
apparent  from  the  -1  g  slopes  in  the  velocity  waveforms. 
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slope.  From  the  figure,  we  can  see  that  the  simulation  somewhat 
over-predicts  the  initial  spall  velocity,  and  therefore  the 
ballistic  period.  This  is  especially  noticeable  near  ground  zero, 
less  so  at  large  range.  The  representative  ballistic  period  near 
ground  zero,  for  the  simulation,  appears  to  be  about  1.7  seconds, 
compared  to  the  observed  period  for  P ILEDRI VER  which  appears  to  be 
about  1.2  seconds.  The  simulation  replicates  very  well  the  decrease 
of  ballistic  period  with  increasing  range,  as  well  as  the  decay  of 
peak  particle  velocity.  Overall,  the  simulation  represents  the  main 
features  of  the  PILEDRIVER  recordings  fairly  well,  and  we  proceed  to 
analyze  the  teleseismic  Rayleigh  wave  from  the  simulation. 

A  one-dimensional  (spherically  symmetric)  explosion  simulation 
was  also  performed,  using  the  same  nonlinear  constitutive  model  as 
was  used  for  the  two-dimensional  simulation.  The  amplitude  spectrum 
of  the  reduced  velocity  potential  obtained  for  this  source  was  shown 
in  Figure  3b. 

We  compute  the  teleseismic  Rayleigh  wave  from  the  near-source 
results  of  the  numerical  simulation,  using  the  procedure  of  Bache, 
Oay  and  Swanger  (1982).  This  procedure  is  analogous  to  the 
surface-integral  formulation  of  Equation  (3),  but  is  exact  rather 
than  asymptotic.  The  surface  of  integration,  I,  is  a  cylinder 
intersecting  the  free  surface  and  enclosing  the  nonlinear  source 
volume,  with  radius  and  depth  of  1209  meters.  As  shown  in  our 
discussion  of  Equation  (3),  any  failure  of  the  near-source  solution 
to  conserve  momentum  will  lead  to  a  long-period  error  in  the  surface 
wave  prediction  by  this  method.  In  principle,  our  finite  difference 
method  conserves  momentum,  but  there  are  several  potential  sources 
of  error.  First,  if  the  velocities  in  the  source  region  (the 
interior  of  I)  are  still  significantly  nonzero  at  the  time  the 
numerical  solution  is  terminated,  the  vertical  traction  integral  on 
Z  will  be  nonzero  at  the  final  time  step,  and  this  will  lead  to  a 
nonzero  spectral  estimate  of  I ^ (m }  at  <n  »  0.  A  second  possibility 
is  that  the  calculation  may  be  run  to  completion,  but  accumulated 
roundoff  errors  might  render  1^(0)  significantly  nonzero.  A  third 


22 


SYSTEMS.  SCIENCE  AND  SOFTWAm 


possioility  is  that  the  initial  stress  field  assumed  at  the  start  of 
the  calculation  may  be  slightly  out  of  equilibrium.  If  such 
disequilibrium  is  small  and  random,  it  may  have  insignificant  effect 
on  1^(0),  but  even  quite  small  systematic  errors  in  the  initial 
stress  field  can  lead  to  substantial  errors  in  1^(0).  This  is 
because  the  effect  of  any  net  forces  on  I  will  be  integrated  over 
the  duration  of  the  simulation. 

To  verify  that  momentum  is  conserved,  to  sufficient  precision 
for  purposes  of  computing  teleseismic  surface  waves,  we  plot  in 
Figure  5  the  total  vertical  momentum  enclosed  by  E,  as  a  function  of 
time.  Also  plotted  is  the  time  integral  of  the  total  vertical  force 
(with  force  taken  relative  to  its  initial  equilibrium  value)  exerted 
on  E  by  the  exterior  continuum.  In  the  absence  of  numerical  errors, 
these  curves  should  coincide.  The  slight  tendency  for  the  momentum 
curve  to  lead  the  impulse  curve  is  an  artifact  of  the  time  series 
processing  to  which  the  stress  histories  on  z  were  subjected  before 
Deing  integrated,  i.e.,  causal  low-pass  filtering  followed  by 
decimation.  Tne  agreement  is  entirely  satisfactory,  and  the 
terminal  value  of  momentum  is  negligible. 

Fundamental-mode  Rayleigh  wave  spectral  amplitudes  were 
computed  at  a  range  of  1000  km  for  the  two-dimensional  explosion 
simulation  using  a  halfspace  with  P  wave  speed  a,  S  wave  speed  6, 
and  density  p  equal  to  6  km/sec,  3.5  km/sec  and  2.7  gm/cm3, 
respectively.  In  the  uppermost  0.2  km,  the  shear  quality  factor 
Qa  was  10,  with  Q  equal  to  300  elsewhere.  The  result  is  shown 
in  Figure  6a,  together  with  the  analogous  result  for  the 
one-dimensional  explosion  simulation.  The  spectral  differences 
reflect  the  influence  of  two-dimensional  source  effects, 
particularly  spall,  on  the  Rayleigh  wave.  As  predicted  by  the 
equivalent-point-force  elastic  model  of  spall,  the  influence  of 
spall  is  most  important  between  about  2  and  5  seconds  period,  and  is 
negligible  at  periods  in  excess  of  10  seconds.  The  maximum 
enhancement  is  about  a  factor  of  2,  and  occurs  at  a  period  of 
approximately  2.5  seconds. 
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igure  5. 


Total  momentum  (positive  up)  in  the  source  volume, 
Vs,  compared  with  the  total  impulse  on  its  boundary, 
E,  for  the  two-dimensional  explosion  simulation.  The 
fact  that  the  curves  are  nearly  coincident  and  approach 
zero  at  late  time  demonstrates  that  the  computation 
properly  conserves  momentum. 
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A  comparison  of  the  Rayleigh  wave  prediction  for  the  simple 
equivalent-point-force  model  of  spall  with  Rayleigh  waves  from  the 
nonlinear  model  is  made  in  Figure  6.  The  solid  curve  in  this  figure 
shows  the  Rayleigh  wave  enhancement  due  to  spall  as  computed  from 
the  nonlinear  simulation;  that  is,  it  plots  the  ratio  of  the  two 
curves  in  Figure  6a.  This  spectral  ratio  is  compared  to  the 
enhancement  factor  I A I  from  the  equivalent-point-force  model  (given 
by  Equation  8),  which  is  shown  as  a  dashed  curve.  Equation  (8)  was 
evaluated  using  the  f( «)  obtained  from  the  one-dimensional  numerical 
simulation.  Both  spall  mass  m$  and  ballistic  period  T$  were 
adjusted  so  that  the  enhancement  factor  agrees  as  well  as  possible 
with  the  spectral  ratio  obtained  from  the  numerical  simulation.  The 
resulting  value  of  T$  -  1.6  seconds  is  in  very  good  agreement  with 
the  ballistic  periods  inferred  from  the  near-source  waveforms  shown 
in  Figure  4.  The  resulting  value  of  m$  «  2.5  x  10^  kg  is  about 
a  factor  of  2.5  greater  than  Viecelli's  mean  empirical  estimate 
(Equation  1).  The  two  curves  are  in  remarkably  good  agreement, 
especially  considering  the  simplicity  of  the  equivalent-point-force 
model.  We  conclude  that  even  when  source-region  rheology  is 
accounted  for  in  detail,  as  in  the  explosion  simulations  analyzed 
here,  spall  makes  no  contribution  to  the  teleseismic  surface  wave 
spectral  amplitude  at  periods  above  10  seconds. 

To  verify  that  spall  has  not  altered  the  long-period  Rayleigh 
waveform,  we  compare  fundamental  mode  synthetic  seismograms  for  the 
one-  and  two-dimensional  explosion  simulations.  The  three-layered 
earth  structure  used  for  the  synthesis  is  shown  in  Table  1.  The 
synthetics,  at  3000  km  range,  and  including  the  long-period  LRSM 
seismometer  response,  are  shown  in  Figure  7.  The  waveforms  for  the 
two  simulations  are  Indistinguishable,  and  peak-to-peak  amplitudes 
agree  within  less  than  half  of  a  percent.  Thus,  for  teleseismic 
surface  waves,  at  periods  greater  than  10  seconds,  we  can  exclude 
the  possibility  of  phase  anomalies,  as  well  as  amplitude  anomalies, 
due  to  spall. 


25 


srsrrMS.  sc/cncc  and  softwakc 


Spectra  Rolio  Rayleigh  Wave  Spectrum 

L  2-0  }  *  I.I-D.I  ( fim  -  sec) 


Frequency  (Hr) 


Figure  6.  (A)  Fundamental-mode  Rayleigh  wave  displacement  spectra 

for  the  one-  and  two-dimensional  explosion  simula¬ 
tions.  (B)  Rayleigh  wave  enhancement  factor  obtained 
analytically  for  the  equivalent-point-force  model, 
compared  to  the  ratio  of  the  one-  and  two-dimensional 
numerical  results  in  (A). 
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TABLE  1 


PROPAGATION  PATH  MODEL  FOR  THE  SYNTHETIC 
SEISMOGRAMS  IN  FIGURE  7 


Layer 

Thickness 

(km) 

a 

(km/ sec) 

0 

(km/sec) 

p 

(gm/cm3) 

QS 

25 

5.5 

3.18 

2.65 

200 

20 

6.4 

3.70 

2.90 

500 

OO 

8.1 

4.68 

3.50 

1000 
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Figure  7.  Synthetic  fundamental-mode  Rayleigh  waves  for  the  one- 
and  two-dimensional  explosion  simulations,  at  a  range 
of  3000  km.  Two-dimensional  effects,  including  spall, 
have  had  no  perceptible  effect  on  the  amplitude  or 
waveform.  Peak  amplitudes  differ  by  less  than  half  a 
percent. 
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VI.  DISCUSSION 


The  main  conclusion  of  this  study  is  that  spall  cannot 
contribute  significantly  to  the  teleseismic  surface  waves  radiated 
by  buried  explosions  at  periods  in  excess  of  about  10  seconds. 
Therefore,  the  spall  phenomenon  cannot  explain  observations  of 
anomalous  surface  wave  amplitude  and  phase. 

It  has  been  shown  that  previous  analytical  studies  which  have 
found  a  potentially  significant  surface  wave  signal  from  spall  are 
based  on  source  descriptions  which  fail  to  conserve  momentum.  This 
momentum  imbalance  exists  for  both  Viecelli 's  (1973)  impulsive-load 
model  and  Harkrider  and  Bache's  (1979)  wave-suppression  source 
model,  and  it  is  directly  responsible  for  the  long-period  surface 
wave  anomalies  predicted  by  those  models. 

The  impulsive-load  model  has  been  modified  to  produce  an 
equivalent-point-force  model  for  spall  which  includes  both  the 
gravitational  and  inertial  forces  associated  with  detachment  and 
closure.  This  model,  which  conserves  momentum,  predicts  surface 
wave  enhancement  at  short  periods,  but  predicts  no  significant 
effect  at  periods  above  about  10  seconds.  This  simple  analytical 
model  thereby  reconciles  the  apparently  conflicting  results  deduced 
from  numerical  simulations  by  Viecelli  (1973)  and  Bache,  Oay  and 
Swanger  (1982),  respecti vely.  Viecelli  observed  enhancement  of  peak 
Rayleigh  wave  amplitude,  because  he  measured  peak  amplitude  very 
near  the  source  (within  the  finite  difference  grid)  at  a  dominant 
period  of  about  3  seconds;  Bache,  et  al_.  found  that  spall  had  no 
effect  on  peak  Rayleigh  wave  amplitude,  because  they  measured  peak 
amplitude  from  synthetic  teleseismic  waveforms  which  have  much 
longer  dominant  periods. 

The  wave-suppression  model  could  be  similarly  modified  to 
conserve  momentum.  This  could  be  accomplished  by  assigning  a 
frequency  spectrum  to  the  suppression  coefficient.  We  have  seen 
that  suppressing  a  fraction  n  of  the  upgoing  wave  introduces  a 
discontinuity  in  vertical  traction  proportional  to  n.  The 
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requirement  imposed  by  conservation  of  momentum  is  that  the 
resulting  vertical  traction  discontinuity  «~z  (Equation  6)  go  to 
zero  at  lo»*j  period  at  least  as  fast  as  w.  Furthermore,  the 
numerical  simulation  indicates  that  virtually  no  momentum  remains 
trapped  in  the  nonlinear  region  longer  than  about  4  seconds 
following  detonation;  therefore,  the  long-period  decay  of  «~z  must 
begin  with  periods  on  the  order  of  4  seconds. 

Under  these  constraints,  only  short-period  upgoing  waves  can 
be  significantly  suppressed  relative  to  downgoing  waves.  With  a 
suppression  coefficient  which  goes  to  zero  at  least  as  fast  as  u  for 
periods  longer  than  several  seconds,  the  wave-suppression  model 
might  provide  a  useful  technique  for  simulating  energy  loss  in  the 
region  between  an  underground  explosion  and  the  surface;  this  region 
is  probably  highly  attenuative,  as  a  result  of  the  tensile  cracking 
associated  with  spall.  The  foregoing  analysis  demonstrates, 

however,  that  such  increased  attenuation  cannot  be  modeled  simply  as 
a  frequency-independent  suppression  of  the  upgoing  wave.  The  same 
reasoning  also  precludes  approximating  near-surface  losses  as  a 
frequency- independent  suppression  of  the  free-surface  reflection 
coefficient.  The  near-surface  environment  may  be  highly  absorptive, 
but  any  seismic  source  must  conserve  momentum  irrespective  of  such 
energy  losses. 

We  note  in  passing  that  Peppin’s  (1977)  model  of  explosion 

cavity  collapse  as  a  downward  force  impulse  is  physically 
unrealizable  for  the  same  reason.  He  argued  that  the  balancing 

upward  impulse  due  to  detachment  of  the  collapse  mass  would  be 
preferentially  dissipated,  since  the  cavity  ceiling  would  consist  of 
highly  shattered  material,  compared  to  the  cavity  bottom.  This 

violation  of  momentum  conservation  might  be  an  acceptable 
approximation  for  modeling  high-frequency  waveforms,  but  it  must  be 
ruled  out  as  a  source  model  for  teleseismic  surface  waves. 
Similarly,  Peppin's  proposed  explosion  source  function,  which 
consists  of  an  upward  impulse  added  to  a  monopole,  is  physically 
unacceptable,  and  will  lead  to  misleading  results  if  applied  to 
teleseismic  surface  waves. 
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The  numerical  simulation  performed  for  this  study  represents 
the  most  complete  multi-dimensional  theoretical  model  of  a  contained 
explosion  of  which  we  are  aware.  It  includes  the  depth  dependence  of 
overburden  pressure,  the  free  surface,  gravity,  depth-varying 
geologic  structure,  and  several  modes  of  material  failure.  The 
calculation  was  time-stepped  well  beyond  spall  closure,  as  verified 
by  the  momentum  curves  in  Figure  5.  Nonetheless,  its  long-period 
Rayleigh  waveform  and  amplitude  are  indistinguishable  from  those  of 
a  spherically  symmetric  explosion,  according  to  Figure  7. 
Furthermore,  the  short-period  Rayleigh  wave  spectral  enhancement 
predicted  by  the  numerical  simulation  is  nearly  identical  in  shape 
to  the  enhancement  predicted  by  the  equivalent-point-force  model 
(Figure  6b).  The  latter  result  is  particularly  surprising  in  view 
of  the  complexity  of  the  numerical  model  and  the  highly  idealized 
nature  of  the  equivalent-force  model.  It  appears  that  the  surface 
wave  prediction  depends  in  a  fairly  simple  way  on  the  total  mass  and 
momentum  of  spalled  material.  Refinements  of  the  nonlinear  material 
models  may  affect  the  amount  of  short-period  enhancement,  but  will 
not  alter  our  conclusion  that  periods  above  10  seconds  are 
uninfluenced  by  spall. 

Clearly,  we  will  have  to  resort  to  phenomena  other  than  spall 
to  explain  surface  wave  anomalies  from  underground  explosions. 
Several  possibility  have  been  proposed  in  the  seismological 
literature,  including  effects  of  deviatoric  prestress  acting  on  the 
source  volume,  (e.g.,  Archambeau,  1972),  explosion-triggered  stress 
release  on  faults  (e.g.,  Ak i  and  Tsai,  1972),  and  explosion-driven 
block  motion  across  rock  joints  (e.g.,  Bache  and  Lambert,  1976; 
Salvado  and  Minster,  1980).  However,  for  none  of  these  mechanisms 
is  the  source  physics  well  understood.  Our  experience  in  studying 
spall  shows  that  multi-dimensional  numerical  simulations,  when 
carefully  analyzed  and  interpreted,  provide  a  useful  tool  for 
studying  complex  source  physics  and  quantifying  its  efre^s  on  the 
seismic  signal. 
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